library(sasld)

# -------------------------------------------------------------------
# Sezione 13.1 - Usare più di una variabile per prevedere una variabile di risposta
# -------------------------------------------------------------------
# Esempio 2
data(hsp)

fit <- lm(Prezzo ~ Dimensione + Stanze, data=hsp)
summary(fit)

fit$coef

model <- fit$model
model[1,]

fit$coef[1] + fit$coef[2]*1679 + fit$coef[3]*3

c(model[1,1],fit$fitted[1],fit$res[1])

a1 <- fit$coef[1] + fit$coef[3]*1
a2 <- fit$coef[1] + fit$coef[3]*2
a3 <- fit$coef[1] + fit$coef[3]*3
c(a1,a2,a3)
b <- fit$coef[2]
plot(c(5000,10000),c(350000,750000),type="n",
xlab="Dimensione",ylab="Prezzo previsto")
abline(a1,b,lwd=2)
abline(a2,b,lty=2,lwd=2)
abline(a3,b,lty=3,lwd=2)
legend(5000,700000,c("1 stanza","2 stanze","3 stanze"),lty=c(1,2,3),lwd=2)

cor(hsp$Prezzo,fit$fitted)

sfit <- summary(fit)
sfit$r.squared




# -------------------------------------------------------------------
# Sezione 13.2 - L'uso della regressione multipla per fare inferenza
# -------------------------------------------------------------------
# Esempio 4
data(college)

fit <- lm(TBW ~ HGT + BF + AGE, data=college)
anova(fit)


ndf <- data.frame(HGT=66,BF=0.18,AGE=20)
ndf

nfv <- predict(fit,new=ndf,interval="prediction")
nfv

nfv <- predict(fit,new=ndf,interval="confidence")
nfv


round(confint(fit),2)




# -------------------------------------------------------------------
# Sezione 13.3 - Analisi dei residui
# -------------------------------------------------------------------

fit <- lm(Prezzo ~ Dimensione + Stanze, data=hsp)

rs <- rstandard(fit)
qqnorm(rs)
abline(0,1)

nok <- which(abs(rs) > 3)
rs[nok]

plot(fit$fitted,rs)




# -------------------------------------------------------------------
# Sezione 13.4 - Regressione multipla con una variabile indipendente qualitativa
# -------------------------------------------------------------------

fit <- lm(Prezzo ~ Dimensione + Condition, data=hsp)
summary(fit)

y <- hsp$Prezzo/1000
x <- hsp$Dimensione*(0.3048^2)
z <- hsp$Condition
fit <- lm(y ~ x + z); summary(fit)

plot(x,y,xlab="Dimensione (metri quadrati)",
ylab="Prezzo (migliaia di dollari)")
ok <- which(z == 1); points(x[ok],y[ok],pch=19)
a <- fit$coef[1]; b <- fit$coef[2]; abline(a,b,lty=3,lwd=2)
a <- fit$coef[1]+fit$coef[3]; b <- fit$coef[2]; abline(a,b,lty=1,lwd=2)
legend(50,850,c("Case in buone condizioni","Case in non buone condizioni"),
lty=c(1,3),lwd=2,bty="n")

fit.a <- lm(y ~ x*z)
summary(fit.a)

a0 <- fit.a$coef[1]
a1 <- fit.a$coef[1]+fit.a$coef[3]
b0 <- fit.a$coef[2]
b1 <- fit.a$coef[2]+fit.a$coef[4]
c(a0,b0)
c(a1,b1)

plot(x,y,xlab="Dimensione (metri quadrati)",
ylab="Prezzo (migliaia di dollari)")
ok <- which(z == 1); points(x[ok],y[ok],pch=19)
abline(a0,b0,lty=3,lwd=2)
abline(a1,b1,lty=1,lwd=2)
legend(0.5,8.5,c("Case in buone condizioni","Case in non buone condizioni"),
lty=c(1,3),lwd=2,bty="n")

anova(fit,fit.a)




# -------------------------------------------------------------------
# Sezione 13.5 - modello logistico
# -------------------------------------------------------------------
# il modello interpolato
x <- seq(0,100,0.1)
plot(c(0,100),c(0,1),type="n",xlab="Reddito",ylab="Probabilità")
lgt <- -3.52 + 0.105*x
p <- exp(lgt)/(1+exp(lgt))
lines(x,p,lwd=2)


# la figura nel testo -----------------------------------------------
x <- seq(0,100,0.1)
plot(c(0,100),c(0,1),type="n",xlab="Reddito",ylab="Probabilità")

lgt <- -3.52 + 0.105*x
p <- exp(lgt)/(1+exp(lgt))
lines(x,p,lwd=2)

lgt <- -3.52 + 0.15*x
p <- exp(lgt)/(1+exp(lgt))
lines(x,p,lwd=2,lty=2)

lgt <- -3.52 + 0.08*x
p <- exp(lgt)/(1+exp(lgt))
lines(x,p,lwd=2,lty=3)

legend(0,1,c("b = 0.08","b = 0.105","b = 0.15"),lty=c(3,1,2),lwd=2)


#
x <- seq(0,100,0.1)
plot(c(0,100),c(0,1),type="n",xlab="Reddito",ylab="Probabilità")

lgt <- -3.52 + 0.105*x
p <- exp(lgt)/(1+exp(lgt))
lines(x,p,lwd=2)

lgt <- -6 + 0.105*x
p <- exp(lgt)/(1+exp(lgt))
lines(x,p,lwd=2,lty=2)

lgt <- -2 + 0.105*x
p <- exp(lgt)/(1+exp(lgt))
lines(x,p,lwd=2,lty=3)

legend(0,1,c("a = -2","a = -3.52","a = -6"),lty=c(3,1,2),lwd=2)
# -------------------------------------------------------------------




# interpolazione
data(cci)

fit <- glm(y ~ income, data=cci, family=binomial(link = "logit"))
summary(fit)


#
si <- c(911, 44, 3,  2)
no <- c(538,456,43,279)
alc <- c(1,1,0,0)
sig <- c(1,0,1,0)
cbind(alc,sig,si,no)

n <- si+no
p <- si/n
cbind(alc,sig,si,no,n,p)

fit <- glm(p ~ alc+sig, family=binomial(link = "logit"), weights=n)
summary(fit)

cbind(alc,sig,si,no,n,p,fit$fitted)


#
fit <- glm(p ~ sig, family=binomial(link = "logit"), weights=n)
summary(fit)

b <- fit$coef
ci <- confint(fit)
exp(cbind(b,ci))

fit <- glm(p ~ alc, family=binomial(link = "logit"), weights=n)
summary(fit)